% Parameters
rho_c = 1050 * 3600;  % Heat capacity [J/(kg*K)]
kappa = 0.564;        % Thermal conductivity [W/(m*K)]
Qz_max = 3.587e-3;    % Max heat source range in z direction
% 假设边界温度为 55°C
% 设置椭圆区域的参数
r_max = 20e-3;   % r方向的最大值（m）
z_max = 20e-3;   % z方向的最大值（m）

% 椭圆的参数：r方向和z方向的半长轴
a = r_max;   % r方向的半长轴
b = z_max;   % z方向的半长轴

% 设置网格的分辨率
n_r = 100;    % r方向的网格点数
n_z = 100;    % z方向的网格点数

% 创建网格
r_vals = linspace(0, r_max, n_r);   % 在r方向上均匀分布的网格点
z_vals = linspace(-z_max, z_max, n_z);   % 在z方向上均匀分布的网格点
[R, Z] = meshgrid(r_vals, z_vals);   % 创建网格的二维矩阵

% 椭圆边界条件：仅保留在椭圆区域内的网格点
ellipse_mask = (R.^2 / a^2 + Z.^2 / b^2) <= 1;  % 计算椭圆区域
R = R(ellipse_mask);  % 仅保留椭圆区域内的r值
Z = Z(ellipse_mask);  % 仅保留椭圆区域内的z值

% 显示椭圆区域内的网格点
figure;
scatter(R, Z);
xlabel('r (m)');
ylabel('z (m)');
title('椭圆区域内的网格点');

% 更新边界温度条件：应用于椭圆边界上的点
for i = 1:length(R)
    if Z(i) == z_max || Z(i) == -z_max || R(i) == r_max
        % 在椭圆边界上设置固定温度
        T(i) = boundary_temp;
    end
end

% Q parameters for forward (z > 0) and backward (z < 0)
alpha_forward = 4.7290;
beta_forward = -3.5090;
C0_forward = 0.3758;
C1_forward = -1.4329;
C2_forward = -0.6797;
C3_forward = 4.8155;

alpha_backward = 4.7290;
beta_backward = -3.5090;
C0_backward = 0.2414;
C1_backward = -0.1028;
C2_backward = -0.1971;
C3_backward = 4.8155;

% Discretization (mesh size)
n_r = 50;              % number of points in r direction
n_z = 50;              % number of points in z direction
n_t = 100;             % number of time steps
dr = r_max / n_r;      % step size in r
dz = 2 * z_max / n_z;  % step size in z
dt = 0.1;              % time step
time_end = 7;          % total time to simulate
% Heat source term
% 重新计算 Q(r, z) 在椭圆区域内的值
Q = zeros(n_r, n_z);  % 初始化 Q 数组

for i = 1:n_r
    for j = 1:n_z
        r = (i-1) * dr;   % 当前的 r
        z = (j-1) * dz - z_max;  % 当前的 z
        
        % 判断是否在椭圆区域内
        if (r^2 / a^2 + z^2 / b^2) <= 1
            if z > 0
                % 正向：z > 0
                Q(i,j) = alpha_forward * exp(beta_forward * r) * ...
                         (C0_forward * z^3 + C1_forward * z^2 + ...
                          C2_forward * z + C3_forward);
            else
                % 反向：z < 0
                Q(i,j) = alpha_backward * exp(beta_backward * r) * ...
                         (C0_backward * z^3 + C1_backward * z^2 + ...
                          C2_backward * z + C3_backward);
            end
        else
            % 超出椭圆区域的点，Q为0
            Q(i,j) = 0;
        end
    end
end


% Initialize temperature matrix
T = ones(n_r, n_z) * 20;  % Initial temperature of 20°C

% Time-stepping loop
for t = 1:n_t
    T_new = T;  % Copy of the current temperature distribution
    
    for i = 2:n_r-1
        for j = 2:n_z-1
            % Finite difference approximation for the heat equation
            d2T_dr2 = (T(i+1,j) - 2*T(i,j) + T(i-1,j)) / dr^2;
            d2T_dz2 = (T(i,j+1) - 2*T(i,j) + T(i,j-1)) / dz^2;
            
            % Time evolution equation
            T_new(i,j) = T(i,j) + dt * (kappa * (d2T_dr2 + d2T_dz2) + Q(i,j) / rho_c);
        end
    end
    
    % Update temperature distribution
    T = T_new;
    
    % Check boundary conditions at z = 5.587e-3, r = 0
    if abs(T(round(n_r/2), round(5.587e-3/dz))) > 55 && abs(T(round(n_r/2), round(5.587e-3/dz))) < 60
        disp('No such solution found: Temperature exceeds 60°C');
        break;
    end
end
% Plot temperature distribution at t = 7s (or other specific time)
figure;
imagesc(linspace(0, r_max, n_r), linspace(-z_max, z_max, n_z), T);
colorbar;
xlabel('r (m)');
ylabel('z (m)');
title('Temperature Distribution at t = 7s');


